Tropical volcanoes synchronize eastern Canada with Northern Hemisphere millennial temperature variability

Although global and Northern Hemisphere temperature reconstructions are coherent with climate model simulations over the last millennium, reconstructed temperatures tend to diverge from simulations at smaller spatial scales. Yet, it remains unclear to what extent these regional peculiarities reflect region-specific internal climate variability or inadequate proxy coverage and quality. Here, we present a high-quality, millennial-long summer temperature reconstruction for northeastern North America, based on maximum latewood density, the most temperature-sensitive tree-ring proxy. Our reconstruction shows that a large majority (31 out of 44) of the coldest extremes can be attributed to explosive volcanic eruptions, with more persistent cooling following large tropical than extratropical events. These forced climate variations synchronize regional summer temperatures with hemispheric reconstructions and simulations at the multidecadal time scale. Our study highlights that tropical volcanism is the major driver of multidecadal temperature variations across spatial scales.

Although global and Northern Hemisphere temperature reconstructions are coherent with climate model simulations over the last millennium, reconstructed temperatures tend to diverge from simulations at smaller spatial scales. Yet, it remains unclear to what extent these regional peculiarities reflect region-specific internal climate variability or inadequate proxy coverage and quality. Here, we present a high-quality, millennial-long summer temperature reconstruction for northeastern North America, based on maximum latewood density, the most temperature-sensitive tree-ring proxy. Our reconstruction shows that a large majority (31 out of 44) of the coldest extremes can be attributed to explosive volcanic eruptions, with more persistent cooling following large tropical than extratropical events. These forced climate variations synchronize regional summer temperatures with hemispheric reconstructions and simulations at the multidecadal time scale. Our study highlights that tropical volcanism is the major driver of multidecadal temperature variations across spatial scales.
Millennial temperatures reconstructed from climate proxies provide crucial historical insights into the temporal and spatial variability of the Earth's climate, as well as benchmarks for quantifying the recent warming and evaluating the realism of climate model simulations [1][2][3][4] . Global and Northern Hemisphere temperature reconstructions now depict high coherence with simulations at the multidecadal time scale, enhancing the predictability of the climate system 2,5 . In contrast, at subcontinental scales, preindustrial temperature variations seem to be region-specific 6 and less spatially coherent in reconstructions than simulations 7,8 . However, it remains unclear to what extent such regional differences reflect high internal variability of the climate system or inadequate proxy coverage and quality [8][9][10] , especially in regions where high-quality proxy data are lacking.
Northeastern North America (hereafter the NENA region) is an archetypal example of this situation where regional uncertainties are large. It is one of the regions of the Northern Hemisphere that lacks millennial-long maximum latewood density (MXD) data, the most sensitive tree-ring proxy for reconstructing summer temperature variability with annual resolution 5,[11][12][13] . Only 12 MXD chronologies have so far allowed reconstruction of summer temperatures back to 1000 CE, with 11 of them clustered in Eurasia (Supplementary Table 1), leading to a well-known gap in MXD-based temperature reconstruction in North America [13][14][15] (Fig. 1a).
In this study, we fill this gap and develop a highly replicated, temperature-sensitive MXD network from four sites along a 900 km latitudinal transect. Based (Fig. 2d), even if this record was developed~450-1400 km away from the sites of our network (Fig. 1a). D-STREC also shares a high fraction of decadal variability (r = 0.50 over the 997-2006 common period, P < 0.001; Fig. 2c) with a previous summer temperature reconstruction of NENA based on independent proxy types (tree-ring width 18 plus δ 13 C and δ 18 O from tree-ring cellulose 19,20 ; Methods), highlighting a well-reconstructed low-frequency domain over the last millennium. The spatial domain of  Table 6).
The "95% CI" is the uncertainty of the Bayesian model, while the "95% CI + data uncertainty" additionally considers the proxy-level uncertainty (Methods). c Comparison of 10-year low-pass filtered (Butterworth filter) D-STREC and 3P-STREC (thick lines) as well as their "95% CI + data uncertainty" (thin lines). D-STREC (significantly positive r with 10-year high-pass CRU temperatures, P < 0.05) covers the northeastern half of North America and has a striking resemblance with the correlation field of the regional MJJA temperature target, even outside the North American continent (Fig. 1b, c). Compared to both 3P-STREC and its tree-ring-width component (Methods), the D-STREC reconstruction skill is substantially improved (Supplementary Table 3), with a much larger spatial domain ( Fig. 1) and better-constrained uncertainties (Fig. 2c), due to the enhanced sensitivity of MXD data to high-frequency temperature variability 5,11,21 (r D-STREC = 0.72 vs. r 3P-STREC = 0.40 with the 10-year high-pass temperature target). Furthermore, D-STREC agrees nicely with the historical record at 1816 CE following the 1815 CE Tambora eruption, although the recovery is one year longer for sites farther north (Supplementary Fig. 3). We speculate that this phenomenon reflects the increasing severity of post-eruption growth stress towards the northern treeline. In contrast, 3P-STREC underestimates this cooling by as much as~2.5°C ( Supplementary Fig. 3a), indicating that volcanic cooling is more reliably recorded by our MXD records. This evidence, in addition to the superior reconstruction skill, implies that D-STREC more precisely reconstructs cold conditions than 3P-STREC, for example, during~1600-1820 CE where the two reconstructions tend to differ (Fig. 2c).
Millennial summer temperature history D-STREC shows that the strongest centennial warming occurred during the twentieth and twenty-first centuries (1917-2016 CE), with a linear temperature rise of 0. 16 22 . The 1590-1930 interval represents an exceptionally cold phase in NENA (−0.69°C), with nine out of the 10 coldest decades and more than half of the 44 coldest summers [≤mean−2 times the standard deviation (SD)] since 772 CE (Supplementary Tables 5, 6).

Impact of volcanic eruptions
Cold decades during the past 1246 years were frequently related to large tropical eruptions. The 1815 CE Tambora eruption resulted in the second coldest summer (1816 CE, −4.53°C with respect to 1905-2006 CE), and the coldest decade (1816-1825 CE, −1.70°C), although this cold period may have been influenced by the 1822 CE Galunggung eruption. Other large tropical eruptions, such as the Huaynaputina (1600 CE), Parker (1640 CE), and Cosiguina (1835 CE) events, were each followed by very cold decades (ranked the 2nd, 4th, and 9th coldest, respectively), although additive effects of multiple events are also possible (Supplementary Table 5). Although the 1453-1462 CE interval is only ranked as the 22nd coldest decade, it is nevertheless the coldest decade between 1000 and 1600 CE, and corresponds to two large tropical eruptions during the 1450 s (1452 and 1457 CE). Conversely, the 1257 CE eruption of Samalas, the most sulfur-rich eruption of the Common Era 23 and associated with a pronounced summer cooling in Europe 24,25 , only produced a weak cooling anomaly of −0.76°C at 1258 CE (with respect to 1905-2006 CE) and a mean temperature of −0.23°C during 1258-1267 CE. This moderate response is in line with the 3P-STREC reconstruction (Fig. 2c) and is probably due to regionspecific volcanic responses 25 and complex atmosphere chemistry [26][27][28] . At the same time, the Samalas event belongs to a sequence of four closely spaced strong tropical eruptions in 1229, 1257, 1275, and 1285 CE (Supplementary Table 7), which may have induced a long-term cooling trend in NENA (−0.06°C per decade during 1200-1300 CE according to D-STREC), at the end of the Medieval Climate Anomaly 18,29 . In fact, at the annual time scale, out of the 44 coldest summers reconstructed by D-STREC, 15 and 16 years can be attributed to tropical and Northern Hemisphere extratropical (NHET) eruption events, respectively ( Fig. 2c; Supplementary Table 6; Methods). The binomialdistribution test indicates a more robust attribution result since 1000 CE than 772 CE (P = 0.021 and 0.051, respectively; Methods). The higher proportion of unattributed cold extremes in the first millennium (5 out of 9) is most likely due to greater uncertainties in ice-core dating of earlier volcanic events 30 . Accordingly, we constrain our subsequent volcano-related analyses to the 1000-2017 period.
Superposed epoch analysis (SEA; Methods) confirms that tropical eruptions induced longer cooling episodes than did NHET eruptions. On average, tropical eruptions caused about a 5-10 years of significant cooling (0.95 significance level) followed by an additional~2-year recovery to the pre-eruption level (Fig. 3a). The cooling peak generally lagged the tropical eruptions by 1 year and corresponded to the forcing peak ( Fig. 3; Supplementary Fig. 4a). In contrast, the cooling effect of NHET eruptions lasted only~1-3 years and was most frequently significant at the year of the eruption ( Fig. 3b; Supplementary Fig. 4b).
Comparisons of SEAs using multiple subsets of volcanic events validate the more persistent cooling following tropical eruptions (Fig. 3a, b; Supplementary Fig. 5), a result not affected by closely spaced eruptions. Analyses on NHET summer temperature reconstructions and simulations ( Supplementary Fig. 6) further confirm the consistency of these results across the Northern Hemisphere. Stratospheric aerosols injected by tropical volcanoes spread poleward with a residence time of 1-3 years 31,32 , while aerosols of NHET eruptions are mainly constrained to 30-90°N with a shorter lifetime 33,34 . Because tropical eruptions influence a larger oceanic domain with high thermal capacity 32,35 , ocean-atmosphere heat exchanges can cool continental summers 36 beyond the direct aerosol forcing, in a more persistent way compared to NHET eruptions (Fig. 3). In contrast to D-STREC, 3P-STREC shows attenuated cooling peaks that lag tropical and NHET eruptions by about 9 and 5 years, respectively ( Supplementary Fig. 7). This behavior most likely results from the strong biological memory of ring-width data 11 , the only high-frequency component of 3P-STREC.  Table 8; Fig. 4c), indicating a high fraction of synchronous multidecadal variability between NENA and hemispheric summer temperatures. The most prominent difference between D-STREC and simulations during the last millennium concerns the impact of the Samalas eruption (Fig. 4b). D-STREC provides complementary proxy evidence that the short-term cooling effect of the Samalas was disproportionately low compared to its amplitude in CMIP5 climate model simulations 8,28,37 .

Region-hemisphere coherence
By averaging 25 members, the CMIP5 multimodel mean largely masks out the unforced internal variability 7,42-44 . Thus, the significant correlations among D-STREC and NHET temperature reconstructions and the multimodel mean simulation imply that externally forced climate variations produce a strong imprint at both regional and hemispheric scales. This is also well supported by strong NENA-NHET coherence in individual CMIP5 simulations (Supplementary Table 8).
Furthermore, D-STREC shows the highest correlations with the volcano-only ensemble among single-forcing simulations of NHET land summer temperatures (Supplementary Fig. 9). This result, along with SEA ( Fig. 3; Supplementary Fig. 5) and correlation analysis after excluding post-eruption years (Fig. 4d, e), points toward volcanism, in particular tropical eruptions, as the main forcing synchronizing NENA and hemispheric summer temperatures. Consequently, in conformity with detection and attribution studies on hemispheric and global millennial temperatures 2,45,46 , our study based on high-quality proxy data highlights the dominant role of tropical volcanism in shaping multidecadal temperature variations across spatial scales. Predictability of multidecadal temperature variability thus remains challenging without any information on future volcanism.

MXD network and chronology development
Our MXD network consists of series from 1668 radii of 1294 black spruce [Picea mariana (Mill.) B.S.P.] trees from four sites across the eastern Canadian boreal forest ( Fig. 1a; Supplementary Table 2), an extratropical region with typically cold/long winters and warm/short summers. L105, L20, and L135 are three newly sampled sites. To ensure data homogeneity 18,47 , we sampled living trees from the lakeshore forests of corresponding lakes where subfossils were collected. New MXD data from these sites were measured from 1-2 radii of each sample using the X-ray densitometric technique (see Wang et al. 48 for details). The dating of millennial chronologies at these three sites was validated using a subfossil wood sample showing a globally coherent cosmogenic 14 C signature at 774 CE 49 . The existing Quex dataset was obtained from the International Tree Ring Data Bank (National Oceanic and Atmospheric Administration, NOAA; Supplementary Table 2) and was corrected for its erroneous location in the metadata and for the cross-dating of one sample (Supplementary Method 1). All MXD measurements were averaged by the tree before subsequent analysis.
After comparing three standardization methods, including the widely applied regional curve standardization 50 , the regionally constrained individual signal-free standardization (RSFi) 51 method was selected to detrend the MXD series at each site (Supplementary Methods 2, 3). The RSFi method efficiently removed nonclimatic signals (e.g., local competition and disturbances) introduced by trees established in different eras, while preserving the long-term variability. In addition, because the MXD series of black spruce trees are known to exhibit heteroscedastic variance 21 , we compared chronologies calculated from ratios and residuals plus power-transformation 52 via a timeefficient linear scaling reconstruction approach 53 . The RSFi ratio chronologies were chosen for the final D-STREC temperature reconstruction due to better overall performance (Supplementary Method 2). Chronology characteristics were assessed by the EPS 16 , rbar, and mean cambial age ( Supplementary Fig. 2).

MXD-based MJJA temperature reconstruction
The D-STREC summer temperature reconstruction was developed using a Bayesian linear regression approach 19 (Supplementary Method 4). Compared to conventional reconstruction methods, the Bayesian approach provides posterior distributions of the climate variable, taking into account individual proxy likelihoods, thus enabling comprehensive uncertainty assessments and improving the reconstruction skill ( Supplementary Fig. 10). The four MXD chronologies showed optimal temperature responses over the MJJA season of the current growing year (Supplementary Fig. 11). Therefore, regional MJJA temperatures were averaged over an area covering our data network (50-60°N, 65-77°W) from the CRU TS 4.03 dataset 17 (Fig. 1a), and the full calibration period was set to 1905-2006 CE, which is the time interval common to the three longest MXD series plus the Two types of confidence intervals were calculated for the D-STREC reconstruction. The first type (referred to as "95% CI") assesses the uncertainties of the Bayesian model and was derived directly from the 2.5th and 97.5th percentiles of the posterior temperatures for each reconstructed year. The second type (referred to as "95% CI + data uncertainty") additionally considers the time-varying uncertainties in the proxy chronologies. We produced 100 chronologies per site with a sampling procedure based on available MXD data points. Assuming that data points of individual trees are normally distributed for each year, the chronologies were built by sampling each year from a normal distribution N~(μ, σ), where μ is the mean of the available data points and σ is the standard error of mean. Supplementary Fig. 1 illustrates the range of ±1.96 × standard error of the mean for the four local MXD chronologies. The 100 chronologies were then included in the Bayesian models for generating alternative D-STREC reconstructions. The final "95% CI + data uncertainty" was derived from the 2.5th and 97.5th percentiles of the mixed posteriors from the 100 runs.

Comparison with the 3P-STREC temperature reconstruction
We compared D-STREC with the earlier 3P-STREC reconstruction developed from published millennial ring width 18 , δ 13 C 19 , and δ 18 O 20 data from the same region near site L20. 3P-STREC and D-STREC are completely independent by proxy types (no proxy in common) and almost independent by sampling sites (one out of 9 sites in common). 3P-STREC comprises one high-frequency component (from ring-width data) and three low-frequency components (ring width, δ 18 O, and δ 13 C data). To allow direct comparison, we recalibrated the 3P-STREC dataset following the same Bayesian procedure and same MJJA temperature target (which correlates significantly with the three individual proxies; Supplementary Fig. 11) as for D-STREC. Similarly, two types of confidence intervals were produced for the 3P-STREC. The ring-width component of 3P-STREC (3P-STREC-TRW), which is a variant of the dataset originally used to develop the first tree-ring-based millennial summer temperature reconstruction in NENA 18 , was additionally compared with D-STREC using a linear scaling approach 53 and the aforementioned regional temperature target.

Historical temperature record
We generated a long MJJA temperature record from a compilation of historical daily temperature observations from the Saint-Lawrence Valley in southern Quebec 54 . The historical record was constructed based on multiple observers in Quebec City, Montreal, and the Ottawa River region, which are~450-800 km away from our closest site, L105 (Fig. 1a). Available maximum and minimum temperature data were averaged to represent the daily mean temperatures, which were later aggregated to monthly data. In order to minimize uncertainties caused by a high frequency of missing values in the early 19th century, a monthly aggregate for each year was retained only if there were less than eight missing daily values in the corresponding month. The valid temperatures for 4 months from May to August were then averaged to yield a seasonal (MJJA) temperature record starting from 1805 CE, with continuous values since 1837 CE (Fig. 2d).

Last-millennium temperature simulations
We used 25 full-forcing and 23 single-forcing last-millennium simulations of monthly near-surface air temperatures. The full-forcing CMIP5 simulations include 16 members from the CESM Last Millennium Ensemble 55 (CESM-LME; including 3 members of isotopeenabled CESM 56 ), and 9 members from the CMIP5 Past1000 41 and corresponding historical 40 experiments. The single-forcing simulations were obtained from the CESM-LME runs singly forced by greenhouse gases, land use, orbital, solar, and volcanic forcing. In addition, the preindustrial 850 control run of the CESM-LME was used as an unforced baseline to evaluate correlations between D-STREC and the single-forcing simulations. Corresponding climate models and experiments are detailed in Supplementary Table 9. In order to allow for direct comparisons among models, all model outputs were interpolated to a T21 resolution (∼5.6°× 5.6°) using the first-order conservative remapping function provided by the Climate Data Operators 57 . Gridded data were then averaged over land between 30°N and 90°N according to area weights to generate a full-forcing and a single-forcing ensemble of simulated NHET land temperatures for the MJJA season. Although full-forcing CMIP5 simulations of Northern Hemisphere summer temperatures tend to underestimate the impact of NHET eruptions (Supplementary Fig. 6), this result probably reflects the fact that CMIP5 volcanic forcing differs from the updated eVolv2k and CMIP6 datasets we used to select volcanic events (Supplementary Tables 7, 9).

Correlation analysis and significance test
We used Pearson's correlation to assess relations among reconstruction, simulation, and climate time series. Because strong autocorrelation reduces effective degrees of freedom of time series and could bias conventional Student's t-tests 58 , we adopted the method of PAGES 2k 2 to test statistical significance for Pearson's r, with a null hypothesis that original time series are unrelated. First, we generated 1000 random red-noise time series for each original series with the same lag-1 autocorrelation coefficient using the "colorednoise" R package 59 . The random series were smoothed, if needed, and then correlated against each other as we did for the original data to form a distribution of 5 × 10 5 correlation coefficients for each pair of comparisons. Finally, these distributions were compared with the true Pearson's r to calculate probabilities (P-values) to test the null hypothesis under a twosided test. Correlation coefficients are considered significant when P is smaller than 0.05.

Attribution of cold extremes to volcanic eruptions
The cold extremes of D-STREC (≤mean−2 SD) were attributed to volcanic eruptions according to three ice-core-based volcanic forcing reconstructions (eVolv2k 30 , IVI2 60 , and ICI 61 ). Locations of eruptions were identified from the corresponding volcanoes provided by the confirmed eruption list (dating uncertainty ≤ 1 year) of the Global Volcanism Program 62 (GVP; tropical: 30°S-30°N, NHET: 30-90°N) or from the ice-core datasets (for unidentified events). Eruptions from Southern Hemisphere extratropics were not considered since they have negligible climatic impacts on the extratropical Northern Hemisphere 63 . Similar to SEA, we updated dates for unidentified eVolv2k events according to Toohey et al. 64 , which are more likely related to true eruption dates (see below). We retained a total of 298 events by removing duplicated eruptions, including those matched by Toohey et al. 30 , from the three eruption datasets. A successful match was identified when a cold extreme corresponded to a retained volcanic event, allowing a maximum 2-year lag prior to the cold year, accounting for time lags of volcanic responses and uncertainties in icecore dating. If several closely spaced eruptions were matched to an extremely cold year, the eruption of the largest magnitude was considered. The attribution result based on the original, date-unadjusted eVovl2k plus IVI2 and ICI datasets is identical to that using the adjusted eVovl2k (Supplementary Table 6). The probability that the observed match frequency differs from a random result was estimated from the binomial distribution. Because of the greater uncertainties in ice-core dating of earlier eruptions 30 , we tested the statistical significance of the attributions for the entire reconstruction and for the period after 1000 CE.

Superposed epoch analysis (SEA)
We used a regular SEA approach provided by the algorithm of Rao et al. 65 to investigate the composite responses to the target volcanic events. First, we constructed a list of tropical and NHET volcanic eruptions (Supplementary Table 7) according to the stratospheric aerosol optical depth at 550 nm (SAOD) estimated from the eVolv2k reconstruction (1000-1900 CE) 30 combined with the Coupled Model Intercomparison Project Phase 6 (CMIP6) volcanic forcing dataset (1901 66 . Events in the first millennium were not considered for SEA due to greater dating uncertainties 30 . The SAOD NHET (areaweighted average over extratropical 30°N to 90°N) ≥ 0.03 (~1/3 of the Pinatubo 1991 CE eruption) was used as criteria by which to select tropical and NHET eruptions that have potentially affected the extratropical Northern Hemisphere. We kept all the events that corresponded to identified eruptions with volcanic explosivity index (VEI) ≥ 4 according to the GVP 62 . The unidentified events were further screened, and were included when also listed in both IVI2 60 and ICI 61 , by permitting a 3-year lag 67 . Four events that fulfilled the criteria (Aira 1471 CE, Serua 1693 CE, Unidentified 1808 CE, and Pinatubo 1991 CE) were discarded because of case-specific reasons (Supplementary Table 7). In total, 24 tropical and 19 NHET events were retained for the SEA.
The key eruption years used for the SEA were re-evaluated to minimize potential uncertainties in the assessment of volcanic effect. For unidentified events, we adopted the years adjusted by the latest application of eVolv2k 64 , accounting for the time lags between eruptions and ice sheet deposition. The adjustments on 10 tropical events led to a more consistent SEA result compared to that using identified tropical eruptions (Supplementary Figs. 5a, 12a). In addition, key years were set 1 year after the eruption for events that occurred in or after August (otherwise assumed in the same year), a strategy adapted from Guevara-Murua et al. 68 . Not considering this lag could introduce biases to SEA because tree ring is a seasonal proxy. For example, an eruption in December cannot affect the tree-ring formation in the same year because black spruce grows in spring-summer 69 . This adjustment was applied to limited events with known eruption months, yet, it could result in a more evident cooling in response to NHET eruptions ( Supplementary Fig. 12b). We performed SEA on multiple time series (reconstructions and simulations) based on the above-constructed key eruption years. We considered the temperature anomalies of 15 posteruption years relative to the 5-year pre-eruption mean. The statistical significance of volcanic cooling was assessed using the block reshuffling method 70 with 10,000 iterations.

Region-hemisphere coherence vs. tropical and NHET eruptions
We designed an additional experiment to compare the impact of tropical and NHET eruptions on the coherence of summer temperatures between NENA and the Northern Hemisphere. We correlated D-STREC with seven NHET temperature series (six reconstructions and the CMIP5 multimodel mean; Fig. 4a, b) after excluding 10 years (the year of eruption included) following each tropical or NHET eruption. This analysis considered the SEA eruption list (key years in Supplementary  Table 7) but excluded the period earlier than 1000 CE and the recent warming trend after 1850 CE. In total, 173 and 143 post-eruption years were excluded from each series to form a nontropical and a non-NHET eruption group, respectively. The retained values in individual temperature series were then chronologically stitched, smoothed using a band-pass filter (20-100 years), and correlated with the similarly processed D-STREC. Although stitching time series is somewhat arbitrary, the resulting correlations help assess the impact of eruptions on the region-hemisphere coherence. We produced 1000 sets of 19 pseudo-tropical and 15 pseudo-NHET eruptions by randomly sampling without replacement between 1000 and 1841 CE, of which the number of post-eruption years to be excluded was the same as corresponding true eruptions. By excluding 10 years following each set of pseudoeruptions, we obtained 1000 correlation coefficients between D-STREC and each NHET series.

Code availability
The code used to analyze data is available from the corresponding author upon request. The code for superposed epoch analysis is available in Rao et al. 65 .